library(tidyverse)
library(censusapi)
library(tigris)
library(sf)
library(mapview)
library(usmap)
library(rjson)
library(censusapi)
library(leaflet)
library(stringdist)
library(measurements)
library(units)
library(pracma)
library(plotly)
library(ggplot2)
library(radiant.data)
Sys.setenv(CENSUS_KEY="b88495f8315f45b2d100dee9ba8f4c489a7371c2")
options(
tigris_class = "sf",
tigris_use_cache = TRUE
)
visits_200309 <- readRDS("P:/SFBI/Restricted Data Library/Safegraph/covid19analysis/weekly-patterns/v2/sj_visits_daily_sum_hourly_03-09-20_05-24-20.rds")
visits_200525 <- readRDS("P:/SFBI/Restricted Data Library/Safegraph/covid19analysis/weekly-patterns/v2/sj_visits_daily_sum_hourly_05-25-20_05-31-20.rds")
visit_data <-
visits_200309 %>%
rbind(visits_200525)
scc_population <-
getCensus(
name = "acs/acs5",
vintage = 2018,
region = "block group:*",
regionin = "state:06+county:085",
vars = c("group(B01003)")
) %>%
mutate(
origin_census_block_group =
paste0(state,county,tract,block_group)
) %>%
rename(
total_pop = "B01003_001E"
) %>%
dplyr::select(total_pop, origin_census_block_group)
sj_blockgroups <- readRDS("sj_blockgroups.rds")
combo_visit_pop <-
visit_data %>%
left_join(
scc_population %>%
filter(origin_census_block_group %in% sj_blockgroups$origin_census_block_group),
by = "origin_census_block_group"
)
date_pull <-
combo_visit_pop %>%
pull(date) %>%
unique()
day_weighted <-
1:length(date_pull) %>%
map_dfr(function(i){
date_filter <-
combo_visit_pop %>%
filter(date == date_pull[i])
single_day_process <-
date_filter %>%
mutate(
multiplied = total_visits_avg * total_pop
) %>%
group_by(date) %>%
summarise(sum_multiplied = sum(multiplied), sum_total_pop = sum(total_pop))
single_day_mean <-
single_day_process %>%
mutate(
weighted_visit_avg = sum_multiplied / sum_total_pop,
mean_percapita = weighted_visit_avg / mean(date_filter$total_pop)
) %>%
dplyr::select(date,weighted_visit_avg,mean_percapita)
single_day_weights <-
date_filter %>%
mutate(
weight_value = total_pop / sum(date_filter$total_pop)
)
single_day_all <-
single_day_mean %>%
mutate(
weighted_visit_sd = weighted.sd(date_filter$total_visits_avg,single_day_weights$weight_value, na.rm = TRUE),
sd_percapita = weighted_visit_sd / mean(date_filter$total_pop)
)
})
block_group_visits <-
combo_visit_pop %>%
mutate(
percapita = total_visits_avg / total_pop
) %>%
dplyr::select(origin_census_block_group,date,percapita)
weekends <-
combo_visit_pop %>%
filter(!duplicated(date)) %>%
arrange(date) %>%
mutate(
x =
case_when(
(date %>% as.numeric()) %% 7 == 1 ~ date + 1,
(date %>% as.numeric()) %% 7 == 4 ~ date - 1,
TRUE ~ date
),
y =
ifelse(
(date %>% as.numeric()) %% 7 %in% c(2,3),
10,
0
)
) %>%
dplyr::select(x,y)
weight_visit_plot <-
plot_ly(
data = day_weighted,
x = ~date
) %>%
add_trace(
x = weekends$x,
y = weekends$y,
type = "scatter",
mode = "lines",
fill = "tozeroy",
fillcolor = "rgba(112,112,112,0.25)",
line = list(color = "transparent"),
name = "Weekends",
showlegend=T
)%>%
add_trace(
y = ~mean_percapita-sd_percapita,
type = "scatter",
mode = "lines",
line = list(color = "transparent"),
name = "City, +/- 1 S.D.",
showlegend=F
) %>%
add_trace(
y = ~mean_percapita+sd_percapita,
type = "scatter",
mode = "lines",
fill = "tonexty",
fillcolor = "rgba(31,164,99,0.25)",
line = list(color = "transparent"),
name = "City, +/- 1 S.D.",
showlegend=T
) %>%
add_trace(
y = ~mean_percapita,
type = "scatter",
mode = "lines",
line = list(color = "rgb(0,0,0)"),
name = "City, Weighted Average",
showlegend=T
) %>%
layout(
margin = list(t = 25),
paper_bgcolor='rgb(255,255,255)',
plot_bgcolor='rgb(229,229,229)',
xaxis = list(
title = "",
gridcolor = 'rgb(255,255,255)',
showgrid = TRUE,
showline = FALSE,
showticklabels = TRUE,
tickcolor = 'rgb(127,127,127)',
ticks = 'outside',
zeroline = FALSE,
fixedrange = T,
range = c(min(combo_visit_pop$date),max(combo_visit_pop$date)+1)
),
yaxis = list(
title = "Mean Visits per Capita",
gridcolor = 'rgb(255,255,255)',
showgrid = TRUE,
showline = FALSE,
showticklabels = TRUE,
tickcolor = 'rgb(127,127,127)',
ticks = 'outside',
zeroline = FALSE,
fixedrange = T,
range = c(0,2)
))
weight_visit_plot